import numpy as np
import matplotlib.pyplot as plt
#根据龙贝格方法计算积分I=1/2pi*(-a到a积分)(1-(1-cosx)/1.2)**c*cosxdx,其中a=arccos(-0.2),c等于1.1或1.5
a=np.arccos(-0.2)
#被积函数为fx=(1-(1-cosx)/1.2)**1.5*cosx
#上下限为-a,a。h=2a/n根据公式T2n=1/2*Tn+h/2*(累加0到n-1)f(x|(k+1/2))
c=0
def f(x):
    fx = (1.0-(1.0-np.cos(x))/1.2)**c*np.cos(x)
    return fx
def leijia(n):
    result = 0
    h = 2*a/n
    for i in range(1,n+1):
        result+=f(-a+(2*i-1)/2*2*a/n)
    return result*h/2
#c=1.5时,为滚珠轴承
c = 1.5
T1 = 1/2*(f(-a)+f(a))
T2 = 1/2*T1+leijia(1)
T4 = 1/2*T2+leijia(2)
T8 = 1/2*T4+leijia(4)
T16 = 1/2*T8+leijia(8)

S1 = 4/3*T2-1/3*T1
S2 = 4/3*T4-1/3*T2
S4 = 4/3*T8-1/3*T4
S8 = 4/3*T16-1/3*T8

C1 = 16/15*S2-1/15*S1
C2 = 16/15*S4-1/15*S2
C4 = 16/15*S8-1/15*S4

R1 = 64/63*C2-1/63*C1
R2 = 64/63*C4-1/63*C2

error = R2-R1

result = R2/(2*np.pi)
print("c=1.5时 龙贝格积分法计算结果为：",result,"误差为：",error)

#c=1.1时,为滚柱轴承
c = 1.1
T1 = 1/2*(f(-a)+f(a))
T2 = 1/2*T1+leijia(1)
T4 = 1/2*T2+leijia(2)
T8 = 1/2*T4+leijia(4)
T16 = 1/2*T8+leijia(8)

S1 = 4/3*T2-1/3*T1
S2 = 4/3*T4-1/3*T2
S4 = 4/3*T8-1/3*T4
S8 = 4/3*T16-1/3*T8

C1 = 16/15*S2-1/15*S1
C2 = 16/15*S4-1/15*S2
C4 = 16/15*S8-1/15*S4

R1 = 64/63*C2-1/63*C1
R2 = 64/63*C4-1/63*C2

error = R2-R1

result = R2/(2*np.pi)
print("c=1.1时 龙贝格积分法计算结果为：",result,"误差为：",error)
#画出函数图像对照
c=1.5
x=np.linspace(-a,a,1000)
y=f(x)*1/(2*np.pi)
plt.plot(x,y)
#传统矩形法求积分验证算法
tryresult = 0
for i in range(1,10001):
    tryresult+=1/(2*np.pi)*f(-a+i/10000*2*a)*2*a/10000
print("c=1.5时 矩形法结果",tryresult)
plt.show()